distrib=function(df,region,level) {
  result=data.frame(matrix(ncol=5,nrow=length(unique(df[[region]]))))
  colnames(result)=c(region,"coef","lb","ub","clusters")
  aa=1
  for (i in unique(df[[region]])) {
    zz=df[df[[region]]==i]
    if (nrow(zz)>500) { 
      if (level=="group") {
        reg=felm(ag2 ~ black.factorvariable | group:year | 0 |group,data=zz) 
      } else if (level == "tract") {
        reg=felm(ag2 ~ black.factorvariable | gtract:year | 0 |group,data=zz) 
      } else if (level == "blockgroup") {
        reg=felm(ag2 ~ black.factorvariable | govbg:year | 0 |group,data=zz) 
      } else { print ("level??") }
      pos=which(substr(names(reg$cse),nchar(names(reg$cse)),nchar(names(reg$cse)))=="B")
      result[aa,1]=i
      result[aa,2]=reg$coefficients[pos,1]
      result[aa,3]=reg$coefficients[pos,1]-1.96*reg$cse[pos]
      result[aa,4]=reg$coefficients[pos,1]+1.96*reg$cse[pos]
      result[aa,5]=nlevels(reg$clustervar[[1]])
      aa=aa+1
    }
  }
  result
}


distrib.county=function(df,region) {
  result=data.frame(matrix(ncol=5,nrow=length(unique(df[[region]]))))
  colnames(result)=c(region,"coef","lb","ub","obs")
  aa=1
  for (i in unique(df[[region]])) {
    zz=df[df[[region]]==i]
    if (nrow(zz)>100 & sum(zz$black.factorvariable=="B")>10 ) { 
      reg=felm(ag2 ~ black.factorvariable | group:year,data=zz) 
      pos=which(substr(names(reg$se),nchar(names(reg$se)),nchar(names(reg$se)))=="B")
      result[aa,1]=i
      result[aa,2]=reg$coefficients[pos,1]
      result[aa,3]=reg$coefficients[pos,1]-1.96*reg$se[pos]  
      result[aa,4]=reg$coefficients[pos,1]+1.96*reg$se[pos]
      result[aa,"obs"]=reg$N
      aa=aa+1
    }
  }
  result
}


distrib.bh=function(df,region,level) {
  result=data.frame(matrix(ncol=5,nrow=length(unique(df[[region]]))))
  colnames(result)=c(region,"coef","lb","ub","clusters")
  aa=1
  for (i in unique(df[[region]])) {
    zz=df[df[[region]]==i]
    if (nrow(zz)>500) { 
      if (level=="group") {
        reg=felm(ag2 ~ blackhisp.factorvariable | group:year | 0 |group,data=zz) 
      } else if (level == "tract") {
        reg=felm(ag2 ~ blackhisp.factorvariable | gtract:year | 0 |group,data=zz) 
      } else if (level == "blockgroup") {
        reg=felm(ag2 ~ blackhisp.factorvariable | govbg:year | 0 |group,data=zz) 
      } else { print ("level??") }
      pos=which(substr(names(reg$cse),nchar(names(reg$cse)),nchar(names(reg$cse)))=="H")
      result[aa,1]=i
      result[aa,2]=reg$coefficients[pos,1]
      result[aa,3]=reg$coefficients[pos,1]-1.96*reg$cse[pos]
      result[aa,4]=reg$coefficients[pos,1]+1.96*reg$cse[pos]
      result[aa,5]=nlevels(reg$clustervar[[1]])
      aa=aa+1
    }
  }
  result
}


distrib.county.bh=function(df,region) {
  result=data.frame(matrix(ncol=5,nrow=length(unique(df[[region]]))))
  colnames(result)=c(region,"coef","lb","ub","obs")
  aa=1
  for (i in unique(df[[region]])) {
    zz=df[df[[region]]==i]
    if (nrow(zz)>100 & sum(zz$blackhisp.factorvariable=="BH")>10 ) {
      reg=felm(ag2 ~ blackhisp.factorvariable | group:year,data=zz) 
      pos=which(substr(names(reg$se),nchar(names(reg$se)),nchar(names(reg$se)))=="H")
      result[aa,1]=i
      result[aa,2]=reg$coefficients[pos,1]
      result[aa,3]=reg$coefficients[pos,1]-1.96*reg$se[pos] 
      result[aa,4]=reg$coefficients[pos,1]+1.96*reg$se[pos]
      result[aa,"obs"]=reg$N
      aa=aa+1
    }
  }
  result
}
